
############# Manuals ##################################
#http://www.bioconductor.org/packages/devel/bioc/html/flowMerge.html

#http://www.bioconductor.org/packages/release/bioc/html/flowQ.html

#http://www.bioconductor.org/packages/release/bioc/html/flowViz.html

#http://www.bioconductor.org/packages/2.10/bioc/vignettes/flowPlots/inst/doc/flowPlots.pdf

#http://www.bioconductor.org/packages/devel/bioc/html/flowTrans.html

#http://www.bioconductor.org/packages/devel/bioc/html/plateCore.html

#http://www.bioconductor.org/help/course-materials/2008/BioC2008/labs/flowCore/flowBasicsSolutions.pdf

################################################################

library(prada);
#pathFiles <- '~/Documents/cytometry/';
pathFiles <- '~/Dropbox/TeamChallenge/cytoData/PAKS/11-11-10-doseresponsetimecourse/';
#cytF <- dir(pathFiles);
#setwd(pathFiles);

cytF <- dir(pathFiles,pattern='*.fcs',recursive=T);

#all_dat <- dat_bu <- read.fcs(cytF);
all_dat <- dat_bu <- readCytoSet(files=cytF,path=pathFiles);


x <- exprs(all_dat[[1]])[,c("FSC-H","SSC-H")];
plot(x, pch = 20, col = densCols(x));


nfit <- fitNorm2(x, scalefac = 2);
plotNorm2(nfit, ellipse = TRUE);

dat_sel <- exprs(all_dat[[1]][nfit$sel,]);


single.well.preprocess <- function(x,scl=2)
{
	dat <- exprs(x);
	fn <- fitNorm2(dat[,c("FSC-H","SSC-H")], scalefac=scl);
	x <- dat[fn$sel,];
	return(x);
}

for (i in 1:length(x))
{
	dat_gated[[i]] <- single.well.preprocess(all_dat[[i]],scl=2);
}



df1 <- list();
dat1 <- list();

for(i in c(1:12))
{
	dat1[[i]] <- subset(drugData[[i]]@data,drugData[[i]]@data[,5]>=mean(drugData[[i]]@data[,5])-3*sd(drugData[[i]]@data[,5]) & drugData[[i]]@data[,5]<=mean(drugData[[i]]@data[,5])+3*sd(drugData[[i]]@data[,5]) & drugData[[i]]@data[,3]>=0 & drugData[[i]]@data[,5]>=0);
	df1[[as.integer(substr(strsplit(strsplit(drugData[[i]]@metadata@filename,'/')[[1]][5],'_')[[1]][3],2,nchar(strsplit(strsplit(drugData[[i]]@metadata@filename,'/')[[1]][5],'_')[[1]][3])))]] <- dat1[[i]][,5]/dat1[[i]][,3];
#	df1 <- list(df1,vec);
}



df2 <- data.frame();
for(i in c(1:12))
{
	vec <- log(df1[[i]],base=2);
	df2[i,1]<-i;
	df2[i,2]<-median(vec);
	df2[i,3]<-mad(vec);
	df2[i,4]<-median(df1[[i]]);
	df2[i,5]<-mad(df1[[i]]);
	df2[i,6]<-mean(df1[[i]]);
	df2[i,7]<-sd(df1[[i]]);
}

tstnls <- nls(V2~(a*V1^n)/(V1^n+k^n)+b,data=df2,start=list(k=10,n=10,a=1,b=.15),trace=T)

k <- coef(tstnls)[1];
n <- coef(tstnls)[2];
a <- coef(tstnls)[3];
b <- coef(tstnls)[4];



png('~/Documents/Hill Fit.png',type='quartz');
plot(df2[,1],df2[,2],pch=16,cex=.7,ylab='log of median fluorescence',xlab='log of est concentrations',main='Hill Equation fit to data');
#curve((7.37*x^10.02)/(x^10.02+7.44^10.02)-4.17,add=T);
curve((a*x^n)/(x^n+k^n)+b,add=T,col='blue');
legend(1,3,c(expression('f(x) = ' (a*x^n)/(x^n+k^n)+b),paste("a = ",round(a,digits=2),"  (",round(confint(tstnls)['a','2.5%'],digits=3),", ",round(confint(tstnls)['a','97.5%'],digits=3),")",sep=''),paste("k = ",round(k,digits=2),"  (",round(confint(tstnls)['k','2.5%'],digits=3),", ",round(confint(tstnls)['k','97.5%'],digits=3),")",sep=''),paste("n = ",round(n,digits=2),"  (",round(confint(tstnls)['n','2.5%'],digits=3),", ",round(confint(tstnls)['n','97.5%'],digits=3),")",sep=''),paste("b = ",round(b,digits=2),"  (",round(confint(tstnls)['b','2.5%'],digits=3),", ",round(confint(tstnls)['b','97.5%'],digits=3),")",sep='')),cex=.85);
dev.off();

plot(df2[,1],df2[,5]/df2[,4],type='l',pch=16,cex=.8,ylab='CV',xlab='log of est concentrations',ylim=c(0.2,1),main='CV Plot');
lines(df2[,1],df2[,7]/df2[,4],ylim=c(0.2,1),col='red');
lines(df2[,1],df2[,7]/df2[,6],ylim=c(0.2,1),col='blue');
legend(10,.9,c("sd/median","sd/mean","mad/median"),col=c('red','blue','black'),lty=1,cex=.7);